%% Simulate Permanent Endowment Shock
% by Jaromir Benes
%
% We simulate the same (endowment) shock in three model objects and compare
% the simulation results: a model with a positive asset position, with a
% negative asset position, and a model with a negative asset position but
% with the sign of the asset position variable flipped. We perform the
% simulations using both of the two methods available in IRIS: "deviations
% from control" and in "full levels".

%% Clear Workspace

clear;
close all;
clc;
irisrequired 20140813;

%% Load Model Objects
%
% Load all three model objects created previously: `m1`, `m2`, and `mflip`.

load read_model.mat m1 m2;
load read_model_flip_sign.mat mflip;

%% Simulate Permanent Endowment Shock
%
% Simulate a 10% permanent endowment shock (shock named `Ey` in the model
% file). Create a steady state database for each model <?sstatedb?>, assign
% the size of the shock <?shock?>, and simulate the shock <?simulate?>.
% These are simulation in full levels. Print the paths for `B` (models `m1`
% and `m2) and `mB` (model `mflip`) on the screen. Note that `B` in model
% `m2` and `mB` in model `mflip` are identical up to the sign.

d1 = sstatedb(m1,1:24); %?sstatedb?
d2 = sstatedb(m2,1:24); %?sstatedb?
dflip = sstatedb(mflip,1:24); %?sstatedb?

d1.Ey(1) = log(1.10); %?shock?
d2.Ey(1) = log(1.10); %?shock?
dflip.Ey(1) = log(1.10); %?shock?

s1 = simulate(m1,d1,1:20,'dbOverlay=',true); %?simulate?
s2 = simulate(m2,d2,1:20,'dbOverlay=',true); %?simulate?
sflip = simulate(mflip,dflip,1:20,'dbOverlay=',true); %?simulate?

[s1.B, s2.B, sflip.mB] %#ok<NOPTS> %?display?

%% Simulate the Same Shock as Deviations from Control
%
% To further understand the logic of log-minus variables, simulate the same
% shock as deviations from control (from steady state). All results for `B`
% and `mB` are positive numbers now: this is consistent with the definition
% of a deviation from steady state for log variables in IRIS: $\hat x_t :=
% x_t / \bar x_t$. A path where a negative log variable becomes larger in
% absolute value is therefore reported as numbers _above_ 1 (see the second
% column in the reported time series <?above1?>), and vice versa.

d1 = zerodb(m1,1:20);
d2 = zerodb(m2,1:20);
dflip = zerodb(mflip,1:20);

d1.Ey(1) = log(1.10);
d2.Ey(1) = log(1.10);
dflip.Ey(1) = log(1.10);

s1 = simulate(m1,d1,1:20,'deviation=',true,'dbOverlay=',true);
s2 = simulate(m2,d2,1:20,'deviation=',true,'dbOverlay=',true);
sflip = simulate(mflip,dflip,1:20,'deviation=',true,'dbOverlay=',true);

[s1.B, s2.B, sflip.mB] %#ok<NOPTS> %?above1?

%% Plot Shock Responses
%
% Plot the results of deviation-from-control simulations for all model
% variables The last graph <?bGraph?> has `B` for models `m1` and `m2`, and
% `mB` for model `mflip`).

sty = struct();
sty.line.lineStyle = {'-','-','--'};
sty.line.marker = {'o','none','none'};
sty.line.lineWidth = {1,1,3};

dbplot(s1 & s2 & sflip,0:20, ...
    {'Y','C','R'}, ...
    'transform=',@(x) 100*(x-1), ...
    'style=',sty);
legend('Model 1','Model 2','Model Flip');

figure();
plot(0:20,[s1.B,s2.B,sflip.mB]); %?bGraph?
grid on;
qstyle(gcf(),sty);
title('B and mB');
legend('Model 1 (B)','Model 2 (B)','Model Flip (mB)');
